import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import simps
from astropy.convolution import convolve, Box1DKernel

E = np.loadtxt('E.txt')

a = np.loadtxt('1p5_50K.txt')
b = np.loadtxt('1p5_200K.txt')
c = np.loadtxt('1p5_600K.txt')

kB = 8.61733e-2 #meV/K
occ1 = 1./(np.exp(E/(kB*50))-1) + 1
occ2 = 1./(np.exp(E/(kB*200))-1) + 1
occ3 = 1./(np.exp(E/(kB*600))-1) + 1

box1 = convolve(a/occ1, Box1DKernel(3))
box2 = convolve(b/occ2, Box1DKernel(3))
box3 = convolve(c/occ3, Box1DKernel(3))

area1 = simps(box1[41:398],E[41:398]) # 41,298,398 for 2, 15, 20 meV
area2 = simps(box2[41:398],E[41:398])
area3 = simps(box3[41:398],E[41:398])

plt.plot(E,box1,color='black',lw=2)  #
plt.plot(E,box2/area2*area1,color='blue',ls=':')
plt.plot(E,box3/area3*area1,color='red',lw=2)

plt.xlim(2,10)
plt.ylim(0,25)
plt.xlabel('Energy (meV)',fontsize=20)
plt.ylabel(r"$\rm{\chi}''$ (a.u.)",fontsize=20)
plt.legend(['50K','200K','600K'],fontsize=20,frameon=False)
plt.title(r'$\mathbf{Q}$=0,1.5,4 (r.l.u.)',fontsize=20)
plt.tick_params(which='both',direction='in',labelsize=20)
plt.minorticks_on()
plt.savefig('Chi1p5_ave_boxcar_norm.pdf',format='pdf',bbox_inches='tight',dpi=600)
plt.show()
